% this code is for Jumps and Post-FOMC Announcement Returns in Currency Markets
% this gives filtered exchange rate (and return) data and jump detection
% results

clear
clc

% 18 countries: 1Australia, 2Brazil, 3Canada, 4Euro, 5Hungary, 6India,
% 7Japan, 8Korea, 9Norway, 10New Zealand, 11Poland, 12Russia, 13Singapore,
% 14South Africa, 15Sweden, 16Switzerland, 17Turkey, 18United Kingdom.

% the sample period is jan 1999 to dec 2018
% base case has 15 min frequency

%% 0. initial set up and filtering

% selecting data frequency: 15 min
npday=96;

% sizing
nc=size(ask,2); % # of currencies (ask: ask quote data matrix)
nt=size(date,1);
ny=max(date(:,8));
nq=max(date(:,7));
nm=max(date(:,6));
nw=max(date(:,5));
nd=max(date(:,4));
% date matrix has numberings for year, quarter, month, week, and day

% mapping different time intervals by using intraday numbering
begendday=nan(nd,2);
begendweek=nan(nw,2);
begendmonth=nan(nm,2);
begendquarter=nan(nq,2);
begendyear=nan(ny,2);

for j=1:nd
    tp=find(date(:,4)==j);

    begendday(j,1)=min(tp);
    begendday(j,2)=max(tp);
end
clear j tp

for j=1:nw
    tp=find(date(:,5)==j);
    begendweek(j,1)=min(tp);
    begendweek(j,2)=max(tp);
end
clear j tp

for j=1:nm
    tp=find(date(:,6)==j);
    begendmonth(j,1)=min(tp);
    begendmonth(j,2)=max(tp);
end
clear j tp

for j=1:nq
    tp=find(date(:,7)==j);
    begendquarter(j,1)=min(tp);
    begendquarter(j,2)=max(tp);
end
clear j tp

for j=1:ny
    tp=find(date(:,8)==j);
    begendyear(j,1)=min(tp);
    begendyear(j,2)=max(tp);
end
clear j tp

% delete inactive quotes (initial check)
temp=ones(nt,nc);

for i=1:nc
    for j=4:nt
        if ask(j,i)==ask(j-1,i) && ask(j,i)==ask(j-2,i) && ask(j,i)==...
                ask(j-3,i)
            temp(j,i)=nan;
        end
        if bid(j,i)==bid(j-1,i) && bid(j,i)==bid(j-2,i) && bid(j,i)==...
                bid(j-3,i)
            temp(j,i)=nan;
        end
    end
    clear j
end
clear i

ask=ask.*temp;
bid=bid.*temp;
clear temp

mid=0.5*(bid+ask);

lr=nan(nt,nc);
for i=1:nc
    for j=2:nt
        lr(j,i)=log(mid(j,i))-log(mid(j-1,i));
    end
end
clear i j

% making indicator matrix for effective returns

ie=nan(nt,nc);

for i=1:nc
    for j=1:nt
        if (lr(j,i)>=0)+(lr(j,i)<0)>0
            ie(j,i)=1;
        end
    end
end
clear i j

% delete 3 consecutively zero returns

for i=1:nc
    for j=3:nt
        if lr(j,i)==0 && lr(j-1,i)==0 && lr(j-2,i)==0
            ie(j,i)=nan;
        end
    end
end
clear i j

% filtering CIP violations
temp=abs(log(frdmid)-log(mid)-0.01*(1/12)*intraintd);
tpmean=nanmean(temp,1);
tpstd=nanstd(temp,1);
tpzscore=zeros(nt,nc);
cipremove=nan(nt,nc);
for i=1:nc
    for j=1:nt
        tpzscore(j,i)=abs((temp(j,i)-tpmean(1,i))/tpstd(1,i));
        if tpzscore(j,i)>6
            cipremove(j,i)=1;
        end
    end
    clear j
end
clear i temp tpmean tpstd

for j=1:nd
    tpcip=cipremove(begendday(j,1):begendday(j,2),:);
    tpsum=nansum(tpcip,1);
    for i=1:nc
        if tpsum(1,i)>30
            ie(begendday(j,1):begendday(j,2),i)=nan;
        end
    end
    clear i tpcip tpsum
end
clear j tpzscore

% winsorize extremes (greater than 9 stdev)
% remove returns with bid-ask bounce back effect
% at first, we need mean and stdev of 15 min return
templr=lr.*ie;
metemp=nanmean(templr);
sttemp=nanstd(templr);
ztemp=zeros(nt,nc);

for i=1:nc
    for j=1:nt
        ztemp(j,i)=abs((templr(j,i)-metemp(1,i))/sttemp(1,i));
        if ztemp(j,i)>9
            ie(j,i)=nan;
        end
    end
    clear j
end
clear i

for i=1:nc
    for j=1:nt-1
        if ztemp(j,i)>7
            if lr(j,i)*lr(j+1,i)<0 && abs(lr(j,i)+lr(j+1,i))<...
                    sttemp(1,i)*0.01
                ie(j,i)=nan;
                ie(j+1,i)=nan;
            end
        end
    end
    clear j
end
clear i

clear metemp sttemp templr ztemp

% remove days with less than 50 observations

for j=1:nd
    temp=ie(begendday(j,1):begendday(j,2),:);
    tpsum=nansum(temp,1);
    for i=1:nc
        if tpsum(1,i)<50
            ie(begendday(j,1):begendday(j,2),i)=nan;
        end
    end
    clear i temp tpsum   
end
clear j

% filtered log return
lre=ie.*lr;

% filtered log excess return
rx=(lr(:,1:nc)+(1/100)*(1/96)*(1/365)*intraintd)*10000;

rxie=nan(nt,nc);    % indicator for effective log excess returns
for i=1:nc
    for j=1:nt
        if ie(j,i)>0
            if (rx(j,i)>=0)+(rx(j,i)<0)>0
                rxie(j,i)=1;
            end
        end
    end
    clear j
end
clear i

rxe=rx.*rxie;

for j=1:nd
    temp=rxie(begendday(j,1):begendday(j,2),:);
    tpsum=nansum(temp,1);
    for i=1:nc
        if tpsum(1,i)<50
            rxie(begendday(j,1):begendday(j,2),i)=nan;
        end
    end
    clear i
end
clear j temp tpsum

rxe=rx.*rxie;

%% 1. constructing market return (dollar factor)

lrm=nan(nt,nc+1);
tempsum=nansum(lre,2);
tempn=nansum(ie,2);

for i=1:nc
    for j=1:nt
        if tempn(j,1)>8
            lrm(j,i)=tempsum(j,1)/tempn(j,1);
            if ie(j,i)>0 && tempn(j,1)>2
                lrm(j,i)=(tempsum(j,1)-lr(j,i))/(tempn(j,1)-1);
            end
        end
    end
end
clear i j
for j=1:nt
    if tempn(j,1)>8
        lrm(j,nc+1)=nanmean(lrm(j,1:nc));
    end
end
clear j tempn

% making indicator for market
lrmie=nan(nt,nc+1);

for i=1:nc+1
    for j=1:nt
        if (lrm(j,i)>=0)+(lrm(j,i)<0)>0
            lrmie(j,i)=1;
        end
    end
    clear j
end
clear i

% for convenience, attach the market return to lre
lre=[lre lrm(:,nc+1)];
lreie=[ie lrmie(:,nc+1)];

% we need day indicator
% intraday and daily indicators for usuable days

iday=ones(nt,nc+1);
for i=1:nc+1
    for j=1:nd
        tempsum=nansum(lreie(begendday(j,1):begendday(j,2),i),1);
        if tempsum<50
            iday(begendday(j,1):begendday(j,2),i)=nan;
        end
    end
end
clear tempsum i j k

% daily series of iday
iieday=nan(nd,nc+1);
for i=1:nc+1
    for j=1:nd
        iieday(j,i)=nansum(iday(begendday(j,1):begendday(j,2),i),1);
    end
end
clear i j

for i=1:nc+1
    for j=1:nd
        if iieday(j,i)>49
            iieday(j,i)=1;
        elseif iieday(j,i)<50
            iieday(j,i)=nan;
        end
    end
end
clear i j

%% 2. Jump Detect & Summary

% jump detection: 95%, 6.65, repeating

% intraday volatility pattern
intvol=zeros(npday,nc+1);
temp=zeros(npday,nc+1);

for i=1:nc+1
    for k=1:npday
        for j=1:nd
            if (lreie(npday*(j-1)+k,i)>=0)+(lreie(npday*(j-1)+k,i)<0)>0
                temp(k,i)=temp(k,i)+lreie(npday*(j-1)+k,i);
                intvol(k,i)=intvol(k,i)...
                    +abs(lre(npday*(j-1)+k,i)*lreie(npday*(j-1)+k,i));
            end
        end
    end
end
intvol=intvol./temp;
clear i j k temp

meanvol=mean(intvol);

adjvol=ones(npday,nc+1);

for i=1:nc+1
    for j=1:npday
        adjvol(j,i)=max(1,intvol(j,i)/meanvol(1,i));
    end
end
clear i j

adjfull=ones(nt,nc+1);

for j=1:nd
    for k=1:npday
        adjfull(npday*(j-1)+k,:)=adjvol(k,:);
    end
end
clear j k

% estimating Bipower variation
bv=nan(nt,nc+1);
bvs=zeros(nt,nc+1);
bvn=zeros(nt,nc+1);
w=156;

for i=1:nc+1
    for j=2:nt
        bv(j,i)=abs(lre(j,i)*lreie(j,i)*lre(j-1,i)*lreie(j-1,i));
    end
end
clear i j

for i=1:nc+1
    for j=w+1:nt
        for k=1:w-2
            if (bv(j-k,i)>=0)+(bv(j-k,i)<0)>0
                bvs(j,i)=bvs(j,i)+bv(j-k,i);
                bvn(j,i)=bvn(j,i)+lreie(j-k,i)*lreie(j-k-1,i);
            end
        end
    end
end
clear i j k

% jump test statistics (L)
lstat=zeros(nt,nc+1);

for i=1:nc+1
    for j=w+1:nt
        if (lre(j,i)*lreie(j,i)>=0)+(lre(j,i)*lreie(j,i)<0)>0 && bvn(j,i)>0
            lstat(j,i)=abs(lre(j,i)*lreie(j,i))...
                /(sqrt(bvs(j,i)/bvn(j,i))*adjfull(j,i));
        end
    end
end
clear i j

lstat(isnan(lstat))=0;

% jump indicator

ji=zeros(nt,nc+1);

for i=1:nc+1
    for j=1:nt
        if lreie(j,i)>0 && bvn(j,i)>10 && lstat(j,i)>6.65
            ji(j,i)=1;
        end
    end
end
clear i j

% repeat the detection procedure once more without jump time points
% remove usuable time points in computing bipower if the time points are
% jumps

ieezero=zeros(nt,nc+1);
for i=1:nc+1
    for j=1:nt
        if lreie(j,i)>0
            ieezero(j,i)=1;
        end
    end
end
clear i j

nie=ieezero-ji;
nji=0*ji;
jp=ji;

jjj=7;  % how many times we want to repeat the jump detection procedure

for q=1:jjj
nie=nie-nji;

nbv=zeros(nt,nc+1);
nbvs=zeros(nt,nc+1);
nbvn=zeros(nt,nc+1);

for i=1:nc+1
    for j=2:nt
        nbv(j,i)=abs(lre(j,i)*nie(j,i)*lre(j-1,i)*nie(j-1,i));
     end
end
clear i j

for i=1:nc+1
    for j=w+1:nt
        for k=1:w-2
            if nie(j-k,i)*nie(j-k-1,i)>0
                nbvs(j,i)=nbvs(j,i)+nbv(j-k,i);
                nbvn(j,i)=nbvn(j,i)+nie(j-k,i)*nie(j-k-1,i);
            end
        end
    end
end
clear i j k

nlstat=zeros(nt,nc+1);

for i=1:nc+1
    for j=w+1:nt
        nlstat(j,i)=abs(lre(j,i)*nie(j,i))/(sqrt(nbvs(j,i)...
            /nbvn(j,i))*adjfull(j,i));
    end
end
clear i j

nlstat(isnan(nlstat))=0;

nji=zeros(nt,nc+1);

for i=1:nc+1
    for j=1:nt
        if nie(j,i)>0 && nbvn(j,i)>10 && nlstat(j,i)>6.65
            nji(j,i)=1;
        end
    end
end
clear i j

jp=jp+nji;
end


jup=zeros(nt,nc+1); % positive jump indicator

for i=1:nc+1
    for j=1:nt
        if jp(j,i)>0 && lre(j,i)>0
            jup(j,i)=1;
        end
    end
end
clear i j

jdn=jp-jup; % negative jump indicator

% to compute moments conveniently,
jpnan=nan(nt,nc+1);
jupnan=nan(nt,nc+1);
jdnnan=nan(nt,nc+1);

for i=1:nc+1
    for j=1:nt
        if jup(j,i)>0
            jupnan(j,i)=1;
            jpnan(j,i)=1;
        end
        if jdn(j,i)>0
            jdnnan(j,i)=1;
            jpnan(j,i)=1;
        end
    end
end
clear i j

% to get jump sizes
js=lre.*jpnan;
jssq=(lre.^2).*jpnan;
jsa=jssq.^(0.5);
jsaup=jsa.*jupnan;
jsadn=jsa.*jdnnan;



















